% SHOLL_TREE   Real sholl analysis.
% (trees package)
% 
% [s, dd, sd, XP, YP, ZP, iD] = sholl_tree (intree, dd, options)
% --------------------------------------------------------------
%
% calculates a sholl analysis counting the number of intersections of the
% tree with concentric circles of increasing diameters. Diameter 0 um is
% 1 intersection by definition but typically 4 points are still output into
% XP...
%
% NOTE ! for loop can be translated into matrix operation easily
%
% Input
% -----
% - intree::integer:index of tree in trees structure or structured tree
% - dd::integer: diameter difference of concentric circles {DEFAULT: 50} OR
%     ::vector: diameter values
% - options::string: {DEFAULT: '-e'}
%     '-s' : show intersections
%     '-3s' : show 3D-intersections
%     '-e' : echo how many double intersections are counted
%     '-o' : count only single intersections
%
% Output
% ------
% - s::vector: sholl intersections
% - dd::vector: diameters
% - sd::vector: double intersections
% - XP,YP,ZP::vectors: coordinates of intersection points
% - iD::vector: index of XP YP ZP in dd.
%
% Example
% -------
% sholl_tree (sample_tree, 20, '-s')
%
% See also dist_tree
% Uses ver_tree dA
%
% intersection points between sphere and line segments directly from
% Paul Bourke 1992
% see http://local.wasp.uwa.edu.au/~pbourke/geometry/sphereline/
%
% the TREES toolbox: edit, visualize and analyze neuronal trees
% Copyright (C) 2009  Hermann Cuntz

function [s, dd, sd, XP, YP, ZP, iD] = sholl_tree (intree, dd, options)

% trees : contains the tree structures in the trees package
global trees

if (nargin < 1)||isempty(intree),
    intree = length (trees); % {DEFAULT tree: last tree in trees cell array} 
end;

ver_tree (intree); % verify that input is a tree structure

if (nargin < 2)||isempty(dd),
    dd = 50; % {DEFAULT diameter difference: every 50 um}
end

if (nargin < 3)||isempty(options),
    options = ''; % {DEFAULT: no option}
end

if numel(dd)==1,
    % if dd is a single value make a vector
    eucl = eucl_tree (intree);
    dd = 0 : dd : (ceil (2 * max(eucl) ./ dd) * dd);
end

[X1 X2 Y1 Y2 Z1 Z2] = cyl_tree (intree); % starting and end points of tree segments
N = length (X1); % number of nodes in tree
% set sphere center to (0, 0, 0)
X3 = X1 (1); Y3 = Y1 (1); Z3 = Z1 (1);
s =  zeros (size (dd));
sd = zeros (size (dd));
XP = []; YP = []; ZP = []; iD = [];
for ward = 1 : length (dd),
    % feed line segments into sphere equation and obtain a quadratic
    % equation of the form au2 + bu + c = 0
    a = (X2 - X1).^2 + (Y2 - Y1).^2 + (Z2 - Z1).^2;
    b = 2 * ((X2 - X1).*(X1 - X3) + (Y2 - Y1).*(Y1 - Y3) + (Z2 - Z1).*(Z1 - Z3));
    c = X3.^2 + Y3.^2 + Z3.^2 + X1.^2 + Y1.^2 + Z1.^2 - ...
        2*(X3 .* X1 + Y3 .* Y1 + Z3 .* Z1) - (dd (ward) / 2)^2;
    squ = b .* b - 4 * a .* c;
    iu = squ >= 0;
    u1 = NaN (N, 1);
    u2 = NaN (N, 1);
    warning ('off','MATLAB:divideByZero');
    u1 (iu) = (-b (iu) + sqrt (squ (iu))) ./ (2 * a (iu));
    u2 (iu) = (-b (iu) - sqrt (squ (iu))) ./ (2 * a (iu));
    warning ('on','MATLAB:divideByZero');
    % when u1 or u2 is in [0, 1] then intersection between segment and
    % sphere exists. When both are in that interval then the segment
    % intersects twice.
    u1 ((u1 < 0) | (u1 > 1)) = NaN;
    u2 ((u2 < 0) | (u2 > 1)) = NaN;
    
    % u1 and u2 are then the solutions on the way from (X1, Y1, Z1) to
    % (X2, Y2, Z2)
    % first add u1 points:
    iu1 = ~isnan (u1);
    Xs1 = X1 (iu1) + u1 (iu1) .* (X2 (iu1) - X1 (iu1));
    Ys1 = Y1 (iu1) + u1 (iu1) .* (Y2 (iu1) - Y1 (iu1));
    Zs1 = Z1 (iu1) + u1 (iu1) .* (Z2 (iu1) - Z1 (iu1));
    XP = [XP; Xs1]; YP = [YP; Ys1]; ZP = [ZP; Zs1];
    iD = [iD; ward*ones(length (find (iu1)), 1)];
    % then u2 points:
    iu2 = ~isnan(u2);
    Xs2 = X1 (iu2) + u2 (iu2) .* (X2 (iu2) - X1 (iu2));
    Ys2 = Y1 (iu2) + u2 (iu2) .* (Y2 (iu2) - Y1 (iu2));
    Zs2 = Z1 (iu2) + u2 (iu2) .* (Z2 (iu2) - Z1 (iu2));
    XP = [XP; Xs2]; YP = [YP; Ys2]; ZP = [ZP; Zs2];
    iD = [iD; ward*ones(length (find (iu2)), 1)];

    s (ward)  = sum (iu1) + sum (iu2);
    sd (ward) = sum (iu1 &~ iu2);
end

s (dd == 0)  = 1;
sd (dd == 0) = 0;

if strfind (options, '-o'),
    s = s - sd;
end

if strfind (options, '-e'),
    if sum (sd) > 0,
        warning ('TREES:wrongcounts',[num2str(sum(sd)) ' segments were counted twice']);
    end
end

if strfind(options,'-s'), % show option
    clf; hold on; shine; HP = plot_tree (intree); set (HP, 'facealpha', 0.5);
    for ward = 1 : length (dd),
        plot (X1 (1) + dd (ward) .* sin (0 : pi / 16 : 2 * pi) / 2,...
            Y1 (1)   + dd (ward) .* cos (0 : pi / 16 : 2 * pi) / 2, 'r-');
    end
    ax = gca; grid on; axis image;
    xlabel ('x [\mum]'); ylabel ('y [\mum]'); 
    ax2 = axes ('Position', get (ax, 'Position'), 'XAxisLocation', 'top', ...
           'YAxisLocation', 'right', 'Color', 'none', 'XColor', 'g', 'YColor', 'g');
    hold on; axis image; xlim (get (ax, 'xlim')); ylim (get (ax, 'ylim'));
    HP = plot (X1 (1) + dd / 2, Y1 (1) + .9 * dd (end) * s ./ (2 * max (s)), 'g-');
    set (HP, 'linewidth', 2);
    clear HP; HP (1) = plot (1, 1, 'r-'); HP (2) = plot (1, 1, 'g-');
    legend (HP, {'sphere diameter', 'intersections count'}, 'location', 'SouthWest');
    set (HP, 'visible', 'off');
    tix = 0 : round (max (s) / 3) : max (s);
    set (ax2, 'ytick', .9 * dd (end) * tix ./ (2 * max (s)), 'yticklabel', num2str (tix'));
    HT = ylabel ('sholl intersections'); set (HT, 'color', [0 1 0]);
end

if strfind (options, '-3s'), % 3D show option
    clf; hold on;
    plot_tree (intree);
    HP = plot3 (XP, YP , ZP, 'r.'); set (HP, 'markersize', 24);
    for ward = 1 : length (dd),
        [XS YS ZS] = sphere (20);
        p = patch (surf2patch ( ...
            X1 (1) + (XS .* dd (ward) / 2), ...
            Y1 (1) + (YS .* dd (ward) / 2), ...
            Z1 (1) + (ZS .* dd (ward) / 2)));
        set (p, 'facecolor', [1 0 0], 'facealpha', 0.01, 'edgecolor', 'none');
        plot (X1 (1) + dd (ward) .* sin (0 : pi / 16 : 2 * pi) / 2, ...
            Y1 (1)   + dd (ward) .* cos (0 : pi / 16 : 2 * pi) / 2, 'r-');
    end
    HP = plot3 (X1 (1) + dd/2, Y1 (1) + zeros (size (dd)), ...
        Z1 (1) + .9 * dd (end) * s ./ (2 * max (s)), 'g.-');
    set (HP, 'linewidth', 4);
    clear HP; HP (1) = plot (1, 1, 'r-'); HP (2) = plot (1, 1, 'r.'); HP (3) = plot (1, 1, 'g.-');
    legend (HP, {'sphere diameter', 'intersections', 'count'});
    set (HP, 'visible', 'off');
    title ('sholl intersections');
    xlabel ('x [\mum]'); ylabel ('y [\mum]'); zlabel ('z [\mum]');
    view (3); grid on; axis image;
end
